AFOSR-TR-  7 8 - 10  2 3 


OPTIMAL  CONTROL  OF  DIFFUSION-REACTION  SYSTEMS 


H.  T.  Banks 
M.  C.  Duban 
J.  P.  Kernevez 


Division  of  Applied  Mathematics 
[Brown  University 
Providence,  R.  I.  02912 


Ddpartement  de  Math^matiques  Appligu^es 
University  de  Technologie  de  CompiSgne 
60206  Compiegne,  France 


May  3,  1978 


Invited  lecture,  International  Conference  on  Applied  Non- 
linear Analysis,  University  of  Texas  at  Arlington, 

April  20-22,  1978. 


This  research  supported  in  part  by  the  National  Science 
Foundation  under  grant  NSF-GP-28931x3 , in  part  by  the 
U.S.  Air  Force  under  contract  AF-AFOSR-76-3092,  in  part  by 
Centre  National  de  la  Recherche  Scientifique  under  grant  38097, 
and  in  part  by  the  University  de  Technologie  de  Compifegne. 


78  06  19  073 


Approved  for  publle  release; 
distribution  unlialted* 


I 

I 


( 


M*  fOMB  OfTioi  07  scxnnrxo 

ntJOM  OF  TRANSUXTTAL  TO  DDC 


MtsumcM  (AFtet 


thlM  ttohnloal  rapart  has  bsen  rsvisw«4  and  la 
•aprovsd  for  public  reloas#  lljf  hr  180-18  (Tbl 
Distribution  is  unllaitsda 


$ 


A,  0*  BLOSI 

Tssbolsal  XafsfMtloQ  Offiaw 


Abstract;  We  discuss  several  formulations  of  optimization 
problems  which  arise  in  a natural  way  in  the  investigation  of 
transport  properties  of  artificial  membranes  and  general  dif- 
fusion-reaction media.  Nonlinear  reaction  velocity  approxima 
tions  dictated  by  reactions  of  interest  to  biochemists  place 
the  problems  in  a class  to  which  one  cannot  apply  the  usual 
computational  techniques  (e.g.  gradient,  conjugate-gradient) 
in  a straightforward  manner.  The  inherent  difficulties,  how 
one  might  circumvent  them,  and  some  of  our  initial  efforts 
towards  development  of  feasible  computational  schemes  are  dis' 
cussed. 
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We  consider  control  problems  governed  by  the  following 
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nonlinear  diffusion-reaction  systems: 


3s  — . 3 s ^ a 


3t 


3x 


1+s+ks^ 


(1) 


3a  3 a 


0<x<l,  0<t<T, 


s(0,t)  = SQ(t)  = 0 


a(0,t)  = aQ(t)  = 0 


(2) 


s(x,0)  = fQ(x)  s(x,T)  = fj^(x) 

a(x,0)  = gQ(x)  a(x,T)  = g^^  (x) . 


(3) 


The  control  of  systems  such  as  (l)-(3)  is  of  importance  in  the 
investigation  of  enzymatically  active  artificial  membranes 
similar  to  those  employed  by  D.  Thomas  and  his  coworkers  in  ex- 
periments at  Universite  de  Technologie  de  Compiegne  (see  [2]  for 
more  details) . In  such  systems  the  variables  s and  a repre- 
sent respectively  normalized  variables  for  substrate  and  activa- 
tor concentrations.  The  nonlinear  reaction  term  in  (1)  is  a 
Michaelis-Menten-Briggs-Haldane  type  (see  Chap.  1 of  [1])  velocity 
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approximation  term  for  a reaction  in  which  one  has  inhibition 
by  excessive  substrate.  The  boundary  conditions  are  those 
appropriate  for  a one  dimensional  diffusion-reaction  medivim  in 
contact  with  a reservoir  (at  x = 0)  and  an  electrode  or  im- 
permeable wall  (at  X = 1)  as  depicted  in  Figure  1. 


i 
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Figure  1 


For  the  nonlinear  system  (l)-(2)  it  can  be  argued  that 
multiple  steady-state  solutions  exist  and  the  initial  and  terminal 


3. 


functions  in  (3)  are  taken  to 
That  is,  i = 0,1  are 


f(0)  = g. 

i 

g(0)  = 


be  distinct  such  steady-states, 
solutions  of 

(4) 

f (1)  = 0 

X 

g^(i)  = 0. 


The  basic  question  we  address  here  is:  Given  the  system 
in  an  initial  steady-state  configuration  (fQ,gQ)  at  time  t = 0, 
how  does  one  use  boundary  controls  to  transfer  the  system 

in  time  0£t<T  to  a second  steady-state  configuration  (fj^,gj^) 
and  do  this  in  an  efficient  manner.  That  is,  there  is  some  cost 
associated  with  adding  (or  deleting)  substrate  and/or  activator 
to  the  system  via  the  boundary  controls  and  one  should  try  to 
minimize  some  measure  of  this  cost  as  the  transfer  from  one  steady- 
state  to  another  is  made.  We  take  as  cost  functional  a measure  of 
the  total  flux  (in  the  L2  sense)  of  s and  a into  the  system 
at  the  boundary  x = 0.  Thus,  we  desire  to  choose  control  func- 
tions BQ,aQ  in  some  control  space  ‘2r(e.g.  L2(0,T))  so  as  to 


}dt 


minimize 


J “ I {|^(0»t)|^  + l|^(C',t)|^ 


(5) 
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subject  to  (l)-(3).  (In  general  the  system  (l)-{3)  need  not  be 
exactly  solvable  for  a given  i = 0,1  (i.e.,  controllability 

questions  arise)  and  one  must  replace  the  above  posed  problem  by 
one  of  transferring  ® terminal  state  close  to 

One  thus  actually  considers  for  both  theoretical  and  computational 
purposes  the  modified  problem  of  minimizing  5 J + 

j I s (x,T)-f^  (x)  I ^ + |a  (x,T)-gj^  (x)  I ^}dx  subject  to  (1)  , (2)  and 

s(x,0)  = f a(x,0)  = gQ(x).) 

The  above  might  appropriately  be  called  a "1-dimensional 
medium"  reaction-diffusion  problem.  An  analogous  "0-dimensional 
medium"  problem  is  of  interest  in  the  event  that  one  has  (i)  reaction 
and  diffusion  separated  within  the  medium  or  (ii)  very  rapid  diffu- 
sion (i.e.,  a well-mixed  medium  for  reaction-diffusion).  The  latter 
assumption  is  valid  in  general  models  for  continuously  stirred  tank 
reactors.  In  the  "0-dimensional  medium"  problem  the  spatial  variable 
is  ignored  and  one  has  as  control  system  (for  s = s(t),  a = a(t)) 


ds  ^ 
dt 

Sq  - s - 

at 

a{aQ-a} 

0 < t_<  T , 

s(0) 

o 

II 

s(T)  = fj^ 

a(0) 

= 9o 

a(T)  = gj^. 

(6) 


(7) 


where  one  still  chooses  the  controls  some  space 

^ of  admissible  policies.  However,  now  the  initial  and  ter-  ) 

minal  states  (frt,g,,),  (fi,gT)  are  constants  which  satisfy  j * 


5. 


where  (s^.a®)  = (s^  (0) .a^ (0) ) , (sQ^aJ)  = (Sq (T) ,aQ (T) ) and 
2 

F(s)  H s/ (1+s+ks  ).  The  cost  functional  is  taken  as 


j = f {|sQ(t)-s(t)  + |aQ(t)-a(t)  |^}dt. 


(9) 


Just  as  in  the  case  of  the  "1-dimensional"  problem,  one  can 
show  that  multiple  steady-states  (i.e.,  solutions  of  (8))  are 
possible  for  the  system  (6) . Also,  one  usually  must  consider  a 
modification  of  the  minimization  problem  since  (6),  (7)  may  not  be 
exactly  solvable  (i.e.,  again  controllability  questions  arise). 

There  are  a number  of  interesting  nontrivial  theoretical 
questions  (controllability,  existence,  uniqueness,  etc.)  associated 
with  the  control  problems  formulated  above  but  we  shall  not  discuss 
those  questions  directly  here.  Our  initial  interest  in  these  pro- 
blems arose  from  an  attempt  to  use  computational  schemes  (i.e., 
software  pac)cages)  in  connection  with  experimental  efforts.  From 
the  descriptions  above  one  mxght  anticipate  this  to  be  a rather 
routine  task  since  the  problems  would  appear  tractable  using 
standard  ideas  from  the  theory  of  boundary  control  of  partial  dif- 
ferential equations  in  the  case  of  the  "1-dimensional"  problems 
(see  [2])  or  those  from  the  theory  of  nonlinear  ordinary  differen- 
tial equation  control  problems  in  the  case  of  the  "0-dimensional" 
problem  (see  (3])  along  with  gradient,  conjugate-gradient  type 
numerical  techniques.  Initial  numerical  experiments  revealed  that 
this  is  not  the  case  and  our  efforts  here  will  be  limited  to  an 
explanation  of  the  difficulties  along  with  suggestions  as  to 
possible  alternative  formulations  which  might  lead  to  problems 
amenable  to  solution  on  the  computer. 
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To  facilitate  discussions  of  the  above-mentioned  dif- 
ficulties it  is  helpful  to  consider  the  quasi-steady-state 
approximation  to  the  "0-dimensional  medium"  problem  (a  similar 
approximation  reveals  the  inherent  difficulties  in  the  "1-dimen- 
sional  medium"  problem) . In  light  of  the  small  transient  times 
found  in  experimental  realizations  of  these  models,  one  can  make 
a plausible  argument  that  the  quasi-steady-state  approximations 
are  reasonable  approximations  to  the  problems  formulated  above. 
We  shall  not  do  that  here  but  turn  instead  to  the  problem  of 
minimizing  J given  in  (9)  subject  to  the  constraint  equations 
(steady-state  approximations  to  (6) ) 


F(s(t))  . 0 

aQ(t)  - a(t)  = 0. 


(10) 


Since  in  this  case  Uq  = a,  we  define  for  convenience  the 
variable  p = aa/(l+a)  and  consider  the  problem  of  minimizing 
J while  transferring  a "state"  X®  = (Sq  (0) , p (0) , s (0) ) to  a 
state  = (Sq (T) , p (T) ,s (T) ) subject  to  the  constraint 

SjjCt)  - s(t)  - p(t)F(s(t))  = 0,  0<t<T.  (11) 

A sketch  of  the  surface  in  (Sq,p,s)  space  described  by  (11) 
is  given  in  Figure  2,  where  one  recognizes  the  well-known  "cusp" 
(catastrophe)  surface  of  Whitney  [5]  and  Thom  [4] . In  Figure  2 


the  folds  in  the  cusp  surface  are  projected  down  into  the  {Sq,p) 
plane  as  the  (infinite)  arcs  containing  CA  and  CB.  We  are 
thus  choosing  control  strategies  (paths  in  the  (Sq,p)  plane) 
which  yield  corresponding  "trajectories"  that  move  on  this  (multi 
valued  in  some  regions)  surface. 

Consider  a problem  which  requires  transfer  of  an  initial 
configuration  to  a terminal  configuration  as  depicted 

in  Figure  2.  Two  possible  distinct  control  strategies 
{ (Sq  (t) , p (t) ) } , 0<t^T,  are  depicted  in  Figure  3. 


Figure  3 
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9. 

It  is  clear  that  two  such  strategies  can  be  made  ar- 
bitrarily close  (using  any  reasonable  measure  of  closeness)  in 
the  (Sq,p)  plane  while  the  corresponding  "trajectories" 

(Sq (t) , p (t) ,s (t) ) , 0<t<T,  lying  on  the  constraint  surface  will 

not  be  close.  The  trajectory  corresponding  to  strategy  1 (see  j 

I 

Fig.  3)  "travels"  along  the  lower  fold  (see  Fig.  2)  while  strategy  j 

j 

2 yields  a trajectory  which  during  the  corresponding  time  "travels" 
along  the  upper  fold  of  the  surface  defined  by  (11) . (The  heavy 
lines  with  arrows  in  Fig.  2 represent  jump  discontinuities  in  s 
for  the  quasi-steady-state  model.  For  the  original  problems,  i.e., 
the  non-quasi-steady-state  models,  these  correspond  to  extremely 
rapid  "motion"  from  trajectories  near  the  lower  surface  to  tra- 
jectories near  the  upper  surface.) 

From  these  considerations  it  is  clear  that  the  trajectories 
for  the  quasi-steady  model  are  not  even  continuous  ^ a function 
of  the  control  strategies  and  hence  it  is  not  surprising  that 
methods  (e.g.  gradient,  conjugate-gradient)  involving  derivatives 
(with  respect  to  controls)  of  the  cost  function  are  troublesome 
when  applied  to  the  problems  governed  by  (l)-(3)  or  (6) -(7). 

Once  one  has  visualized  the  problems  in  this  heuristic  but 
informative  way,  it  is  apparent  that  the  difficulties  are  a result 
of  the  particular  nonlinear  reaction  velocity  approximation  found 
in  (1)  and  subsequent  associated  versions  of  this  system  equation  em- 
ployed above.  The  models  entail  a region  F (for  (6)  and  (11) 
with  transfer  from  X®  to  as  shown  in  Figs.  2,3  this  region 

is  depicted  in  Fig.  4)  in  "control"  space  in  which  one  must  choose 
control  strategies  with  extreme  care. 


Figure  4 


In  carrying  out  laboratory  experiments,  this  region  is 
observed  to  be  one  in  which  the  system  is  highly  unstable.  Thus 
from  both  a theoretical  and  practical  viewpoint,  additional  con- 
straints on  operation  of  the  system  in  this  region  are  desirable. 
Careful  formulation  with  additional  constraints  can  lead  to 
tractable  problems.  We  illustrate  this  first  with  a sketch  of 
how  one  might  formulate  such  a control  problem  for  a discretized 
version  of  the  quasi-steady  approximation  to  the  "0-dimensional 
medium"  problem. 


! 
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Considering  the  controls  to  be  piecewise  constant 

on  [0,T],  one  can  reformulate  the  quasi-steady  problem  as  a 
multi-stage  discrete  control  problem  with  "controls" 

{Sq  (ti) ,aQ  (ti) } , i = l,...,k,  constrained  to  lie  outside  r 
(see  Fig.  4)  with  "states"  {s(t^)}  given  implicitly  by 


ao(ti) 


=0<‘i'  - F(s(tj))  = 0. 


The  payoff  is  then  taken  as 


J 


I tSQ(tQ) 
1=1 


- s{t^)}^At^. 


The  most  natural  formulation  along  these  lines  leads  to  immediate 

difficulties  with  regard  to  necessary  conditions  (multiplier 

rules  or  maximum  principles  are  not  readily  available  for  discrete 

control  problems  with  implicit  state  equations) . However,  one 

can  reformulate  this  slightly  as  a constrained  "state"  and  "control" 

problem  so  that  necessary  conditions  are  easily  obtained.  If  one 

2 1 

identifies  ^q'^q  "states"  and  defines  a mapping  A:R  -*■  R 
by  = A(Xj^,X2)  where  x^  is  a solution  (appropriately  chosen 
when  multiple  solutions  exist)  to 


- *3  - "rFi^  ''‘*3’  ■ “ 

and  introduces  "controls"  v^^,  (with  suitable  constraints), 

the  problem  becomes  one  of  minimizing 


( 
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J = {sQ(t^)  - A(Sq (t^) ,aQ (t^) ) }^At^ 


subject  to  state  equations 


and  constraints 


= =0<Vl>  ''i 


* "i 


'®0  ^ ° 


defined  to  prohibit  values  of  SQ,aQ  in  the  region  r.  In  this 
formulation  one  can  obtain  necessary  conditions  (to  use  as  a 
basis  for  computational  schemes)  via  application  of  the  operator 
theoretic  optimization  framework  with  abstract  multiplier  rule 
developed  by  Meustadt  (e.g.,  see  Chap.  7 of  [3]). 

The  above  formulation  essentially  involves  the  assumption 
that  "changes"  (or  more  precisely  "rates  of  changes")  in 
are  the  controls.  This  can  be  viewed  as  a special  case  of  a re- 
formulation for  the  continuous  version  problems.  Consider  the 
full  "0-dimensional  medium"  problem  and  adjoin  to  the  state 
equations  (6)  additional  equations 
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with  control  constraints  lv|  < Mj^,|w|  < M2.  The  "states"  for 
the  problem  are  then  taken  as  SQ,aQ,s,a  with  "controls"  given  by 
v,w.  In  addition  to  the  natural  control  restraints,  one  imposes 
mixed  state-control  (so-called  "phase-control"  constraints)  in- 
equality constraints  which  restrict  the  choices  of  v and  w in 
the  event  one  is  in  the  region  F in  (Sq,p)  space  (see  Fig.  4). 
These  constraints  are  defined  so  that  one  rules  out  control  policies 
that  yield  paths  in  the  (Sq,p)  plane  that  travel  along  the 
"singular"  arc  containing  CA  (see  Figs.  2,3,4).  That  is,  one 
rules  out  via  constraints  on  v,w  policies  such  as  those  depicted 
in  Fig.  3.  Hence,  while  one  does  not  prohibit  crossing  of  the 
r region,  one  restricts  carefully  the  types  of  trajectories  one 
allows  while  passing  through  this  region.  The  resulting  con- 
straints will  thus  be  joint  in  the  "states"  SQ,aQ  and  the  "controls" 
v,w. 

Finally,  we  point  out  that  one  can  also  use  these  reformu- 
lation ideas  to  take  a linear  programming  approach  (function  minimi- 
zation problems  with  linear  inequality  constraints)  to  these  pro- 
blems as  opposed  to  the  optimal  control  approach  (multiplier  rule 
for  constrained  control  problems)  sketched  above.  For  example 
we  illustrate  briefly  with  the  quasi-steady  approximation  to  the 
"0-dimensional  medium"  problem. 

We  approximate  the  arc  containing  CA  above  and  below  by 
straight  lines  with  slope  m.  With  this  fixed  slope  m,  we  con- 
struct  a family  of  parallel  and  equidistant  lines  {p  "'^0^^i^i=0 
(Bq,p)  plane  as  depicted  in  Figure  5.  The  admissible 


in  the 


(Sq{T) ,p{T) ) 


Figure  5 


"states"  are  then  required  to  lie  on  these  lines.  (The  construc- 
tion is  made  so  that  the  arc  CA  lies  between  two  of  these  lines, 
(Sq(0),  p(0))  lies  on  p = mSg  + b^,  and  (Sq(T),  p (T) ) lies 
on  p = mSp  + bj^.)  A "trajectory"  will  then  consist  of  a sequence 
of  points  (SQ(t^),  p(t^),s(t^))  satisfying  (11)  with  (sQ(t^), 
p(t£))  belonging  to  the  line  p = ms^  + b^.  One  can  take  as  con- 
trol policies  the  collection  of  sequences 


U t f p 2r  ^2*  ’ ' * > * ^]c“l ^ 


with  tp  = 0,  tj^  = T and  corresponding  "state”  equations 

So(ti)  = (p^-b^)/m 

SjjCt^)  - s(t^)  - pj^F(s(t^))  = 0 . 

In  addition  to  making  appropriate  modifications  to  the  payoff 
J,  one  constrains  the  analogues  of  equations  (12),  i.e., 


(13) 


One  must  also  add  positivity  constraints  for  Sq,p  given  by 

p^  > max{0,b^}.  (14) 

By  defining  suitable  coefficient  matrices  E and  D,  one  can 
write  the  constraints  (13),  (14)  as 

Eu"^  > D.  (15) 


^i-^i-1 


Pi-Pi_l 


^i-^i-1 


The  problem  then  becomes  one  of  minimizing  a function  J subject 
to  the  linear  inequality  constraints  (15)  and  standard  computa- 


16. 


tional  techniques  le.g.,  descent  methods  such  as  the  Davldon- 
Fletcher-Powell  schemes )are  applicable. 

A more  detailed  discussion  of  theoretical  aspects  of 
the  above  different  formulations  and  approximations  along  with 
our  numerical  findings  will  be  presented  in  a forthcoming  manu- 
script. 
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